/*=================================================================
 *
 * Chebftd.c	.MEX file corresponding to chebftd.M
 * Computes Chebyshev polynomial fit to a set of data...
 *
 * The calling syntax is:
 *
 *		c = chebftd(a, b, n, t, d)
 *
 *  You may also want to look at the corresponding M-code, which is commented out
 * in expfit.m.
 *
 * This is a MEX-file for MATLAB.  
 * Paul B. Manis, Ph.D.
 * 9/13/02 
 *=================================================================*/
/* $Revision: 1.0 $ */
#include <math.h>
#include "mex.h"

/* Input Arguments */

#define	A_IN	prhs[0]
#define	B_IN	prhs[1]
#define N_IN	prhs[2]
#define T_IN	prhs[3]
#define D_IN	prhs[4]


/* Output Arguments */

#define	C_OUT	plhs[0]

#if !defined(MAX)
#define	MAX(A, B)	((A) > (B) ? (A) : (B))
#endif

#if !defined(MIN)
#define	MIN(A, B)	((A) < (B) ? (A) : (B))
#endif

#define PI 3.14159265



static void chebftd2(
		   double	*c,
		   double	a,
 		   double	b,
		   int	n,
		   double	*t,
		   double	*d,
		   int		nd
		   )
{
	double bma, bpa, y, fac, sum, *f;
	double inc;
	int pos;
	int k, j;

	f = (double *) mxCalloc(n, sizeof(double));
	
//	 mexPrintf("\nEntering Chebyshev\n");
	
	bma = 0.5*(b - a);
	bpa = 0.5*(b + a);
	inc = t[1]-t[0];
	for (k=0; k < n; k++) {
		y = cos(PI*(k+0.5) / n);
		pos = (int)(0.5+((y*bma+bpa)/inc));
		if(pos < 0) {
			pos = 0;
		}
		if(pos > (nd-2)) {
			pos = nd-2;
		}
		*(f+k) = d[pos+1];
//	mexPrintf("pos = %d   d(pos+1) = %f   f(%d): %f\n", pos, d[pos+1], k, *(f+k));
	}
	fac = 2.0 / (float) n;
//	mexPrintf("fac: %f bma=%f bpa=%f inc=%f n=%d\n", fac, bma, bpa, inc, n);
	for (j=0; j<n; j++) {
		sum=0.0;
		for (k=0; k<n; k++) {
			sum += f[k]*cos(PI*j*(k+0.5)/n);
		}
		*(c+j)=fac*sum;
//	mexPrintf("j = %d, sum = %f, c[j] = %f\n", j+1, sum, *(c+j));
	}
	mxFree(f);
}

/* Mex function wrapper for chebftd.c */
void mexFunction( int nlhs, mxArray *plhs[], 
		  int nrhs, const mxArray*prhs[] )
     
{ 
    double *c; 
    double a, b, *t, *d; 
	int n, nd;
    unsigned int i,j,k,m; 
    
    /* Check for proper number of arguments */
    
    if (nrhs != 5) { 
	mexErrMsgTxt("chebftd.dll: Five input arguments required."); 
    } else if (nlhs > 1) {
	mexErrMsgTxt("chebftd.dll: Too many output arguments."); 
    } 
     
    /* Create a matrix for the return argument */ 
    n = (int)mxGetScalar(N_IN);
   C_OUT = mxCreateDoubleMatrix(1, n, mxREAL); 
    
    /* Assign pointers to the various parameters */ 
    c = mxGetPr(C_OUT);
    
    t = mxGetPr(T_IN); 
    a = mxGetScalar(A_IN);
    b = mxGetScalar(B_IN);
    d = mxGetPr(D_IN);
	nd = mxGetN(D_IN);
    /* Do the actual computations in a subroutine */
    chebftd2(c, a, b, n, t, d, nd); 
//    mexPrintf("c[1] = %f\n", c[1]);

    return;
    
}


